Flat bands and Wigner crystallization in the honeycomb optical lattice 
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We study the ground states of cold atoms in the tight-binding bands built from p-orbitals on a 
two dimensional honeycomb optical lattice. The band structure includes two completely flat bands. 
Exact many-body ground states with on-site repulsion can be found at low particle densities, for 
both fermions and bosons. We find crystalline order at n — | with a v3 x \/3 structure breaking 
a number of discrete lattice symmetries. In fermionic systems, if the repulsion is strong enough, we 
find the bonding strength becomes dimerized at n = \. Experimental signatures of crystalline order 
can be detected through the noise correlations in time of flight experiments. 
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Optical lattices have opened up a new venue in which 
to study strongly correlated systems, with the possibil- 
ity of precisely controlled interactions. For example, 
the superfluid-Mott insulator transition [l[ for bosons 
has been observed. Many interesting states with novel 
magnetic and supcrfiuid properties in optical lattices 
have also been proposed by using high-spin bosons and 
fermions [2|. In addition, optical lattices provide the 
opportunity for studying another important aspect of 
strongly correlated systems, orbital physics, which plays 
an important role in metal- insulator transitions, super- 
conductivity, and colossal magneto-resistance [|[. 

Cold atom systems with orbital degrees of freedom ex- 
hibit new features which are not usually realized in solid 
state systems. Recently, the properties of the bosons in 
the first excited p-orbital bands have attracted a great 
deal of theoretical attention @, H, H, S B 111 i including 
predictions of a nematic superfluid state [5j, antiferro- 
magnetic ordering of orbital angular momentum (OAM) 
00], a striped phase of OAM in the triangular lattice Q, 
and a bond algebraic liquid phase [9(. Experiments car- 
ried out by Browaeys et al. [l(| and Kohl et al. [ll[ have 
demonstrated the population of higher orbital bands in 
both bosons and fermions. Recently, Foiling et al. has 
showed the existence of both the Mott-insulating and su- 
perfluid states of the p-band bosons [12] • 

The physics of graphene captured by the 2D honey- 
comb lattice has generated tremendous interest as a re- 
alization of Dirac fermions (l3j . In graphene, the active 
bands near the Fermi energy are "tt" type, composed of 
the p z orbitals directly normal to the graphene plane. 
The other two p-orbitals (j> x ,y) lie in-plane, and exhibit 
both orbital degeneracy and spatial anisotropy. How- 
ever, they hybridize with the s-band and the resulting 
CT-bonding band is completely filled. In contrast, in opti- 
cal lattices, p XjV bands are well separated from the s-band 
with negligible hybridization, giving rise to the possibil- 
ity of interesting new physics. The honeycomb optical 
lattice has already been realized quite some time ago, by 
using three coplanar laser beams fl4| . 



In this article, we study the p^y-orbital physics in the 
2D honeycomb lattice. We find the lowest energy band is 
completely flat over the entire Brillouin zone (BZ). When 
the flat band is partially filled, the effects of interactions 
are entirely non-perturbative. For sufficient low densi- 
ties the ground states are "Wigner" crystals (we slightly 
abuse this nomenclature to apply it generally to both 
bosons and fermions). We obtain the exact many-body 
plaquette Wigner crystal state at filling (n) = |. For 
fermionic systems, we obtain additional crystalline or- 
dered states at higher commensurate fillings. The noise 
correlation in the time of flight image should detect the 
order in these states. 

We first discuss the single-particle spectrum. The op- 
tical potential on each site is approximated by a 3D 
anisotropic harmonic well with frequencies uj z 3> oj x — 
Lu y = Lu xy , and thus the p z orbital is well separated in en- 
ergy from the p x j,-orbitals. Since the honeycomb lattice 
is bipartite, we denote by A and B its two sublattices. We 
define three unit vectors ei,2 = ± 2 ^ + \e v , e 3 = — e y 
pointing from each A site to its three B neighbors. The 
projection of the p-orbitals along the three 6^1^,3 direc- 
tions are defined as pi = (p x e x + P y e y ) • e-i- Only two of 
them are linearly independent. The kinetic energy reads 

H o = h ipkiPf+aea + h - c -} - m W 

where the cr-bonding term iy describes the hopping be- 
tween p-orbitals on neighboring sites parallel to the 
bond direction, a is the nearest neighbor distance, and 
np = n? x + npy is the total particle number in both 
p x and p y orbitals at the site r. t\i is positive due to 
the odd parity of the p-orbitals and is set to 1 below. 
Eq. [1] neglects the 7r-bonding t± which describes the 
hopping between p-orbitals perpendicular to the bond 
direction. Typically tn/t± ^> 1 due to the high spatial 
anisotropy of the p-orbitals. We have numerically con- 
firmed this for the optical potential realized in experi- 
ment V(r) = Vq X«=i~3 C0S (Pi ' r) where p, = p e; and 




FIG. 1: Dispersion of the two-lowest p x ^-orbital bands i?i,2- 
The band Ei is completely flat, while E2 exhibits Dirac points 
at Ki^2 = (zfc , 0). The other two bands are symmetric 
with respect to E = 0. 

Po = & Namely, t ± » 0.02t|| and t|| « 0.24Er at 
Vo/i?fl = 5 [E R is the recoil energy). 

The band structure contains both flat bands and Dirac 
cones. Each unit cell consists of two sites, each of which 
contains two orbitals p x . y , resulting in four bands. The 
BZ takes the shape of a regular hexagon with edge length 
47r/(3\/3a). The dispersion of all the bands is symmetric 
with respect to the zero energy because of the bipar- 
tite nature of the lattice. The b and structure con sists of 
Ei A {k) = =F|,S 2i3 (fc) = Ti\/3 + 2]r\ cosfc-6, where 

h = a (e 2 - e 3 ), b 2 = a (e 3 - ei) and b 3 = a (ei - e 2 ). 
Wie show only the two lowest bands E\ t2 in Fig. [TJ 
which touch at the Brillouin zone center. The bottom 
and top bands turn out to be completely flat. The two 
middle bands are dispersive with a width determined by 
f II , which are the same as in graphene with two non- 
equivalent Dirac points located at — (±^^-,0). We 
will not repeat the extensively studied Dirac cone physics 
here, but rather focus on the new features brought by the 
flat bands instead. 

The degeneracy of all momentum eigenstates in the flat 
bands allows us to take any superposition of these states, 
and in particular, to construct eigenstates that are local- 
ized. As depicted in Fig. [2]A, for each hexagon plaquette 
R, there exists one such eigenstate for the bottom band 

6 

= £( : ' : {^''jl'j-r ^n0 !llj:/ } (2) 
i=i 

where j is the site index and 0j — (j — The or- 

bital configuration on each site is perpendicular to the 
links external to the hexagonal loop. This orbital orien- 
tation, together with destructive interference of hopping 
amplitudes, prevent the particle from "leaking" off of the 
plaquette. The states IV'r) are all linearly-independent 
apart from one constraint I^Pr) = under periodic 
boundary conditions. The Bloch wave states are con- 

structed as |V*>a = dr Efc e <£J W(* ^ (°' ^' with 
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FIG. 2: A) The Wannier-like localized eigenstate for the 
lowest band. The orbital configuration at each site is ori- 
ented along a direction tangential to the closed loop on which 
the particle is delocalized. B) The configuration of the close- 
packed Wigner-crystal state for both bosons and fermions at 
71 = 1/6. Each thickened plaquette has the same configura- 
tion as in A). 

the normalization factor Nk = |(3 — ^ cos/c • hi). The 

doubly degenerate eigenstate at k = (0, 0) can not be 
constructed from the above plaquette states. They are 

l^iS=(0,0))l.2 = J2feA \Px(y),r) ~ EfeB \Px(y).r)- 

Because of the orbital degeneracy, the onsite interac- 
tion for spinless fermions remains Hubbard-like as 

Hint — U^^TL^ x Tl^y. (3) 

r 

In order to enhance U between neutral atoms, we may 
use the 53 Cr atom which has a large magnetic moment 
of 6/is (Bohr magneton), and polarize the spin with an 
external magnetic field. The length scale of p xl/ -orbitals 
l x . y = y^h/mujx^y is typically one order smaller than a. 
For example, we estimate that l x< y/a> ~ 0.2 at Vq/Er = 5. 
Assuming strong confinement in the z-axis l z <C l x ,yi the 
vector linking two atoms in p x and p y orbits almost lies in 
the plane. U can be adjusted from repulsive to attractive 
by tuning the polarization direction from perpendicular 
to parallel to the xy-plane. We will only discuss repulsive 
U below. We estimate that U can easily reach the order 
of Er which is much larger than tu . The off-site interac- 
tions are small and decay with distance as 1 /r 3 , and thus 
are neglected. For example, the nearest neighbour inter- 
action is at the order of (^-) 3 U w 10~ 2 U. The on-site 
interaction for p-band bosons is given in Ref . [f| [§[ . 

When the flat band is partially filled the effect of in- 
teractions is non-perturbative. At sufficiently low par- 
ticle density n < 1/6, both the interactions and the ki- 
netic energy can be minimized simultaneously. This can 
be realized if the individual particles localize into the 
plaquette states depicted in Fig. [2] A. These are exact 
eigenstates of the flat band, and therefore minimize ki- 
netic energy. Any arrangement of the plaquette states 
where they do not overlap will cost no interaction energy. 
These localized plaquette states therefore constitute ex- 
act many-body ground states. The maximum density 
(close packed) arrangement of non-overlapping plaque- 
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FIG. 3: The phase boundary of the incompressible plaquette 
Wigner crystal state of spinless fermions at (n) = ~. 



ttes has filling n = 1/6, and has the structure depicted 
in Fig. [2] B. The \/3 x \/3 structure, breaks the lat- 
tice translation symmetry and is three-fold degenerate. 
At zero temperature, the particle density jumps from 
to 1/6 as the chemical potential is increased through 
/.t = —3/2. Similar flat band behavior has been studied 
in the magnon spectra of a number of frustrated magnets 
in a large magnetic field [l5[ near full polarization of the 
magnet. The degenerate states for the Kagome antifcrro- 
magnet are in exact one to one correspondence to those of 
our model. Thus the two models have identical thermo- 
dynamics (if fluctuations are restricted to these states, as 
appropriate for k B T <C tn, U). This is described by the 
classical hard hexagon model[lj|, which exhibits a sec- 
ond order thermal phase transition in the 3-state Potts 
model universality class, breaking translational symme- 
try, when the fugacity of the hexagon z c = (11 + V / 5)/2. 
In this state, the atoms do not touch each other, thus 
particle statistics do not play any role. The Wigner crys- 
tal is also expected to appear for bosons in this optical 
lattice at filling n = ^. In contrast, Wigner crystal is not 
possible in graphene systems [l6j . 

When the filling n > 1/6, exact solutions are no longer 
available for the interacting Hamiltonian. For simplicity, 
in this regime we shall only discuss the case of spinless 
fermions, within a mean field treatment. We decouple 
Eq. [3] in both the direct and exchange channels as 



H m f,int — — ^2 { n r,x (»V) — (?V,3) 



+ («r,3>] -PpxPfy - *(«f,2) — /l.C.}, (4) 

and solve it self-consistently. Here rai,a,3 are the 
pseudo-spin operators defined as n? 3 = \{p\ x Pf x — 

pl y Pfy),nf,l = \{P ?x Pry + h.C.),Jl^ 2 = ^{p^Pfy ~ fl.C.). 

The Tift t i, 7if?3 operators are time-reversal invariant, and 
describe the preferential occupation of a "dumbbell- 
shaped" real p-orbital orientation; np 2 is the orbital an- 
gular momentum, and is time-reversal odd. At the mean 
field level (np^) is zero. In order to obtain the plaquette 
order in Fig. [5J we take the six sites around a plaquette 
as one enlarged unit cell in the mean field calculation. 
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FIG. 4: The filling (n) vs. the chemical potential fi for spinless 
fermions for weak A) and strong B) interactions. Due to 
particle-hole symmetry, only the part with from the band 
bottom — §t|| to (7/2 is shown. Only one plateau appears 
in A) at n = ~, while a series of plateaus appear in B) at 
n= 1/6, 1/3, 1/2,2/3,5/6, 1. 



The numerically determined phase boundary between 
the 1/6 state and a compressible phase at higher density 
is shown in Fig. [3] The charge gap grows roughly linearly 
with U in the weak interaction regime, and saturates at a 
value comparable to iy in the strong interaction regime. 
This can be understood as follows: in the weak interac- 
tion case, we choose a plaquette R which is adjacent to 
three occupied plaquettes ^1,2,3 as depicted in Fig. [2j B, 
and put an extra particle in it. The cost of the repulsion 
is ^r. In the strong coupling case, we put the particle 

into an excited state of the occupied plaquette R± while 
fixing the orbital configuration on each site. Because 
fermions are spinless, the cost of energy comes from the 
kinetic part with the value of |iy. Thus the charge gap 
A < min(i?7, |ty) which agrees with the numeric result. 

The curves of the filling n vs. /i in both the weak and 
strong coupling regimes are depicted in Fig. [4] In both 
cases, the plaquette Wigner state appears as a plateau 
with (n) = g. In the weak coupling regime (U/t\\ = 1), 
as [i passes the charge gap, (n) increases quickly corre- 
sponding to filling up other states in the flat band. Due 
to the background crystal ordering, those states are no 
longer degenerate, but develop weak dispersion. A signif- 
icant reduction in density of states (DOS) occurs at the 
approximately commensurate filling of (n) » ^, but it is 
not a strict plateau. At (n) > 1/2, after the flat band 
is completely filled, the interaction effects are no longer 
important. Near half-filling, the DOS vanishes where the 
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FIG. 5: The dimerized state of spinless fermions with filling 
(n) — |. Each thickened (red) bond corresponds a dimer 
containing one particle. 

physics is dominated by the Dirac cones. 

The physics changes dramatically in the strong cou- 
pling regime. We looked specifically at (U/t = 10), where 
we found a series of plateaus appear at commensurate fill- 
ings (n) = g(i = 1 ~ 6). The large charge gap at (n) = 1 
is of the order of U. The other plateaus with (n) < 1 are 
the novel features the p-orbital model presents. In ad- 
dition to the plaquette bond order at (n) — i discussed 
above, rich structures including dimer and trimer bond 
orders appear at other plateaus with (n) < 1. We will 
present these patterns in detail in a future publication. 
Here a dimer refers to a superposition of the two states 
of two sites where one is occupied, while the other is 
empty. The dimerized state at (n) = h is illustrated in 
Fig. O Each dimer is represented as a thickened bond 
where the orbital configuration is along the bond. We 
check that the bonding energy for the thickened bond is 
approximately 0.95t|| at U/t\\ = 10 which is about one 
order larger than that of other bonds. 

Next we discuss the effect of adding a small 7r-hopping 
t± term. Then the lowest energy band acquires weak 
dispersion, with a small band width of t±. This associ- 
ated kinetic energy cost for the n = | plaquette state 
is of order t± per particle. A stability condition of this 



state can therefore be roughly estimated as ^ > t±. We 
have checked numerically that the ^-plateau survives at 
U > U with the value of t± = O.ltn. Other plateaus 
appearing in the strong coupling regime are not sensitive 
to small t±. 

The plaquette Wigner crystal phase in Fig . [2]B should 
manifest itself in the noise correlation for the time of 
flight (TOF) signals. In the presence of the plaquette or- 
der, the reciprocal lattice vector for the enlarged unit cell 
becomes G x = (^,0) and G 2 = (^, The cor- 

relation function is defined as Ct(r, , r f ) = {n{r)n{f l )) t — 
(n(f)}t(n(r'))t- After a spatial averaging and normaliza- 
tion, we find 

C t (d) = [ dr Ct ^+^ ? -i\ K ± y 6 (k - G)(5) 

where '+' ('— ') is for bosons (fermions) respectively, k — 
md/{ht), and G = mG\ + nGi with m,n integers. 

There are numerous directions open for further explo- 
ration. Some interesting variations on the model are to 
consider spin-fui fermions, for which there is a possibil- 
ity of ferromagnetism as a result of the flat bands, or 
attractive interactions, in which pairing and the BCS- 
BEC crossover in the flat band might prove interesting. 
Most intriguing is the possibility of exotic incompressible 
states analogous to the Laughlin liquid in FQHE. These 
cannot be captured within the mean-field approximation 
used here for n > 1/6. If one could devise appropri- 
ate variational liquid states projected into the flat band, 
these could be compared energetically with the Wigner 
crystals found here. Given the richness and surprises en- 
countered in the FQHE, flat band physics in optical lat- 
tices appears rife with possibility. The comparison with 
the fermion condensation where the flat band arises due 
to strong interactions will also be interesting [171 ]. 
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